Quantum Monte Carlo calculation of the equation of state of neutron matter 
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We calculate the equation of state of neutron matter at zero temperature by means of the auxiliary 
field diffusion Monte Carlo method (AFDMC) combined with a fixed-phase approximation. The 
calculation of the energy is carried out by simulating up to 114 neutrons in a periodic box. Special 
attention was made to reduce finite size effects at the energy evaluation by adding to the interaction 
the effect due to the truncation of the simulation box, and by performing several simulations using 
different number of neutrons. The finite size effects due to the kinetic energy were also checked by 
employing the twist-averaged boundary conditions. We considered a realistic nuclear Hamiltonian 
containing modern two- and three-body interactions of the Argonne and Urbana family. The 
equation of state can be used to compare and to calibrate other many-body calculations and to 
predict properties of neutron stars. 



I. INTRODUCTION 

The equation of state of nuclear matter and its prop- 
ertiesplays a central role in the modeling of neutron 
stars The density in the star ranges from a small 
fraction of, up to several times, nuclear saturation den- 
sity, 0.16 fm^"^, which is found in the center of heavy 
nuclei. At such extreme conditions no phenomenologi- 
cal data determined from experiments are available, and 
because the matter inside the neutron stars is closer to 
neutron matter than symmetric nuclear matter, heavy- 
ion collision experiments do not substantially constrain 
the equation of state A realistic calculation of the 
equation of state of neutron matter is then particularly 
challenging in both many-body nuclear physics and as- 
trophysics. 

The equation of state of neutron matter can in prin- 
ciple be computed in the framework of many-body the- 
ories using a bare interaction. A common alternative 
is represented by effective Skyrme forces. However, the 
resulting equation of state strongly depends on the pa- 
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rameters of the effective interaction used, even in the low 
density regime At present, there are a wide range and 
type of Skyrme interactions. However their non-realistic 
character impairs their ability to reliably calculate the 
properties of neutron stars|4i]. More accurate many-body 
techniques are then needed to perform predictive calcu- 
lations. 

A microscopic calculation of neutron matter starting 
from a non-relativistic nucleon-nucleon and three- nucleon 
interaction is both challenging and of great relevance. 
Variational techniques based on correlated basis func- 
tions are good candidates to solve for the ground-state of 
neutron matter. The operatorial structure of the nuclear 
Hamiltonian and the strong correlations induced by the 
high density make these techniques hard to use. The en- 
ergy evaluation using the correlated basis function theory 
is usually performed by solving the Fermi Hyper Netted 
Chain (FHNC) equationsQ neglecting many elementary 
diagrams. In addition, the operatorial structure of the 
Hamiltonian leads to additional approximations, like the 
Single Operator Chain (SOC) approximation, due to the 
non-commutativity of the terms entering in the varia- 
tional wave function. Therefore the resulting equation 
of state contains, in principle, uncontrolled approxima- 
tions which may be partially corrected by computing the 
energy exactly up to a few first orders in the cluster ex- 
pansion 0. 

Despite the progress of the last several years in the 
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determination of sophisticated two- and three-nucleon 
interactions, large discrepancies among different calcu- 
lations of nuclear and neutron matter are still present. 
Quantum Monte Carlo techniques based on projection 
can be very accurate for calculating the ground state and 
low lying excited states of nuclei. In particular Green's 
Function Monte Carlo (GFMC) was employed to fit the 
three-nucleon interaction form in nuclei up to ^ = 80, 
and then used to test the nuclear Hamiltonian up to the 
^^C ground-state Q. At present the huge number of re- 
quired numeric operations limits the applicability of this 
method to only about 14 neutrons[9]. 

The Auxiliary Field Diffusion Monte Carlo 
(AFDMC[l3]) combined with a fixed-phase approx- 
imation was employed to predict properties of nuclei in 
very good agreement with the GFMC[ll|, and stressed 
important limitations of other many-body theories used 
in nuclear matter calculations 

In this work we present an accurate evaluation of the 
equation of state of neutron matter using a realistic two- 
and three-nucleon interactions (the Argonne v'g and Ar- 
gonne vis combined with Urbana-IX[13]). The computed 
equation of state can be used as a benchmark for other 
many-body techniques. 

The plan of the paper is the following: in the next sec- 
tion we describe the structure of the nuclear Hamiltonian 
we used; in Sec. IIIII we will briefly review the AFDMC 
method and explain the fixed-phase approximation. The 
results will be presented in Sec. llVl and some conclusions 
will be given in the last section. 



II. HAMILTONIAN 

Properties of a generic nuclear system can be studied 
starting from the non-rclativistic Hamiltonian 
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which includes the kinetic energy operator, a two-nucleon 
interaction Wy, and a three-nucleon interaction Vi 



ijk- 



The nucleon-nucleon interactions are usually depen- 
dent on the relative spin and isospin state of the nucleons 
and therefore written as a sum of several operators. The 
coefficients and radial functions that multiply each oper- 
ator are adjusted by fitting experimental scattering data, 
and the type and number of these operators depends on 
the interaction. A large amount of empirical informa- 
tion about the nucleon-nucleon scattering problem has 
been accumulated. In 1993, the Nijmegen group ana- 
lyzed all nucleon-nucleon scattering data below 350 MeV 
published in physics journals between 1955 and 1992 14\. 
Nucleon-nucleon interaction models that fit the Nijmegen 
database with a /Ndata ^1 are called "modern" which 
include the Nijmegen modelsfl^ (Nijm93, Nijm I, Nijm 
II and Reid-93), the Argonne 

models [I6|,ll2| and the CD- 

Bonn[18]. However all of these interactions, when used 



alone, underestimate the triton binding energy, suggest- 
ing that the contribution of a three-nucleon force is es- 
sential to reproduce the physics of nuclei. 

The most sophisticated Argonne interaction is the Ar- 
gonne wistlfi] potential, written as a sum of 18 operators. 
However we often consider another interaction, the Ar- 
gonne WgflT'l that is a simplified version of Argonne wis; 
it contains only 8 operators, and the prime symbol indi- 
cates that such potential is not just a simple truncation 
of Argonne iiig, but a reprojection, which preserves the 
isoscalar part in all S and P partial waves as well as in 
the ^Di wave and its coupling to ^Si. The Argonne Wg 
is a bit more attractive than Argonne wig in light nuclei 
by about 0.5 MeV per nucleon[13j|. but its contribution 
is very similar to Argonne wig in neutron drops, where 
the difference is about 0.06 MeV per neutron[7]. 

The Argonne potential between two nucleons i and j 
is written in the coordinate space as a sum of operators 



(2) 



0=1 



where n is the number of operators which depends on 
the potential, Vp{r) are radial functions, and r^j is the 
inter-particle distance. 

The eight operators included in Argonne Wg give the 
largest contributions to the nucleon-nucleon interaction. 
The first six of them come from the one-pion exchange 
between nucleons, while the last two terms depend on the 
velocity of nucleons and give the spin-orbit contribution. 
These eight operators are 

Of~^'® = (IjCTi • crj,Sij,Lij ■ Sij) X (1,T,; • Tj) , (3) 

where Sij is the tensor operator 

Sij = 3((T, • r.y)(crj • f,j) - (Ti ■ (Tj , (4) 

Lij is the relative angular momentum of couple ij 

L^J^Y^^n-rJ)xiV^-Vj), (5) 
and Sij is the total spin of the pair 



1 



(6) 



with both Lij and Sij divided by h to make them unit- 
less. 

In modern interactions these eight operators are the 
standard ones required to fit S and P wave scattering 
data in both triplet and singlet isospin states. 

The three-nucleon interaction contribution is mainly 
attributed to the possible A intermediate states that an 
excited nucleon could assume after and before exchanging 
a pion with other nucleons. This process can be written 
as an effective three-nucleon interaction, and its param- 
eters are fit to light nuclei 0, [l§| and eventually to prop- 
erties of nuclear matter, such as the enipirical equilib- 
rium density and the energy at saturation[20|. The three- 
nucleon interaction must accompany the two-nucleon in- 
teraction and the total Hamiltonian studied. 
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The Urbana IX three nucleon interaction used in our 
calculation has the following form 

V,,u = + Ffl . (7) 
The Fujita-Miyazawa term[2l| is spin-isospin dependent: 



where 



Y{x) 
T{x) 



-[Xij,X.jk\[Ti ■ Tj,Tj ■ Tk] 



i + - + ^)n^)eT(r) 



(8) 
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The phenomenological Vr part is 
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eye 
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For neutrons, the commutator terms in Eq. [5] are zero, 
and each of the anticommutator terms has only spin op- 
erators for two of the three neutrons. 

The A2ti term of Urbana-IX was originally fitted, along 
with the Argonne wig parameters, to reproduce the triton 
and alpha particle binding energy, while the Uq strength 
was adjusted to obtain the empirical equilibrium den- 
sity of nuclear matter^2|- However, while the ground 
state of light nuclei can be exactly solved with few-body 
techniques, the determination of the equation of state 
of symmetric nuclear matter can be evaluated only using 
many-body techniques that contain uncontrolled approx- 
imations. 



III. METHOD 
A. Diffusion Monte Carlo 

The Auxiliary Field Diffusion Monte Carlo method is 
an extension of the usual Diffusion Monte Carlo to deal 
with Hamiltonians that are spin-isospin dependent. The 
diffusion Monte Carlo method [23l. [24 |. projects out the 
ground state properties by starting from a trial wave 
function not orthogonal to the true ground state. 

Consider a generic trial wave function ipx expanded 
over a set {0n} of eigenstates of the Hamiltonian: 



VtCR) =^(i?,0) = ^c„(/.„(R) 



(11) 



The propagation in imaginary time r is given by 

V'lR, r)=J2 c„e~(^~^«'^<^„(R) , (12) 



where Eq is a normalization factor, and R represent the 
spatial coordinates of the system. In the limit t —^ oo, 
iPCRjt) approaches the lowest eigenstate (po with the 
same symmetry as ip. The evolution can be done by 
solving the integral equation 



V'(R,t) = y"G(R,R',T)V(R',0)dR', 



(13) 



where the wave function is described with a set of 
configurations called walkers as following: 



and 



(R|^) = ^(R)-£(R|Rfe)(Rfe| 

fc=i 



(R|Rfc) = ^(R-Rfe). 



(14) 



(15) 



The kernel G'(R, R', r) is the Green's function of the sys- 
tem, and can be expressed as the matrix clement 



G(R,R',r) = (R|e-(^-^o)nR'> 



(16) 



By considering a generic Hamiltonian and the Trotter 
decomposition, the form of the Green's function in the 
small imaginary time-step limit At —^ is 



G(R, R', At) « e — 2 ' 



^Go(R,R',At), (17) 



where Go is the Green's function of the noninteracting 
system 



(18) 



and the factor due to the interaction plus the trial eigen- 
value Et is the normalization of the Green's function 
factor computed over the time interval r: 

' V{R) + V(R') „^ 

(19) 



-Bt ) Ar 



The integral equation [T3] can be solved in a Monte 
Carlo way. At each time-step all walkers are moved with 
the diffusion term of the free Green's function Go, so that 
for each walker a new set R' of spatial coordinates are 
generated according to 



R' = R + 77 , 



(20) 



where R is the old configuration, and 77 is a vector of 
random numbers with probability density Gq. 

The normalization of Eq. [191 translated into a weight 
of the walker, is sampled using the branching technique 
in which w gives the probability of a configuration to 
multiply at the next step. Computationally, this is im- 
plemented by weighting estimators according to w, and 
generating from each single walker a number of replicas 



n — [w + S] 



(21) 



where ^ G [0; 1] is a random number and [x] means integer 
part of X. 

The infinite imaginary-time limit is reached by iterat- 
ing this process for a sufficient total time r = nAr. 
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B. Auxiliary Field Diffusion Monte Carlo 

In the case of nuclear Hamiltonians the potential con- 
tains quadratic spin and isospin and tensor operators, 
so the many-body wave function cannot be written as a 
product of single particle spin-isospin states. 

For instance, let us consider the generic quadratic spin 
operator CTj • <jj where the cr are Pauli's matrices oper- 
ating on particles. It is possible to write 



2P° 



(22) 



where P^^ interchanges two spins, and this means that 
the wave function of each spin-pair must contain both 
components in the triplet and singlet spin-state [25l [2^. 
By considering all possible nucleon pairs in the systems, 
the number of possible spin-states grows exponentially 
with the number of nucleons. 

Thus, in order to perform a diffusion Monte Carlo cal- 
culation with standard nuclear Hamiltonians, it is nec- 
essary to sum over all the possible single particle spin- 
isospin states of the system to build the trial wave func- 
tion used for propagation. This is the standard approach 
in GFMC calculations for nuclear systems. 

The idea of AFDMC is to rewrite the Green's function 
in order to change the quadratic dependence on spin and 
isospin operators to a linear dependence by using the 
Hubbard-Stratonovich transformation. 

For neutrons Ti ■ Tj — 1, so that the isoscalar-spin 
operators of the Hamiltonian can be recast in a more 
convenient form 



V - 5]5]^^p(r„)0(P)(z, j)l^s/ + VsD 

ia,jf3 



(23) 



where Latin indices label nucleons, Greek indices label 
Cartesian components, and 



Vsi = ^ [vi{rij) + V2{i 



(24) 



is the spin-isospin independent part of the interaction. 
The 3A by 3A matrix v4'^°'^ contains the interaction be- 
tween nucleons of other terms; 

Mnj) + v,in,)] (Sfr^rf^ - 5^p) . (25) 

The matrix A is zero along the diagonal (when i — j),in 
order to avoid self interaction, and is real and symmetric, 
with real eigenvalues and eigenvectors given by 



j75 



(26) 



The matrix A^"'^ has n = 1...2>A eigenvalues and eigenvec- 
tors. We can then define a new set of operators written 
in terms of eigenvectors of the matrix A: 



The spin-dependent part of Eq. [23] becomes 
2 



and the corresponding propagator is then 
At first order in At this is equivalent to 



n- 



(27) 



(28) 



(29) 



(30) 



Each factor can be linearized with respect to the oper- 
ators O by using the Hubbard-Stratonovich transforma- 
tion 



1 



'2tt J- 



dxe 
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(31) 



The Green's function then becomes: 



G(R,R',Ar) = ( 



2Trh^AT 

3A 



-Vs7(R.)At 



^\ V2^ J 



(32) 



The newly introduced variables a;„ , called auxiliary fields^ 
are sampled in order to evaluate the integral of Eq. [321 
For a spin state given by a product of single particle spin 
functions, the linearized Green's function has the effect 
of changing the spin state by independently rotating the 
spin of each single nucleon. 

The sampling of auxiliary fields to perform the integral 
of Eq. [31] eventually gives the same effect as the propaga- 
tor with quadratic spin operators acting on a trial wave 
function containing all the possible good spin states. The 
effect of the Hubbard-Stratonovich is then to reduce the 
dependence of the number of operations needed to eval- 
uate the trial wave function from exponential to linear 
in the number of nucleons. The price to pay is the addi- 
tional computational cost due to the diagonalization of A 
matrices and the sampling of the integral over auxiliary 
fields. 

Sampling of auxiliary fields can be achieved in several 
ways. The most intuitive one, in the spirit of Monte Carlo 
sampling is to consider the Gaussian in the integral Eq. 
[5T] as a probability distribution. The sampled values are 
then used to determine the action of the operators on the 
spin part of the wave function. This is done exactly as 
in the diffusion process. It is possible to use other tech- 
niques to evaluate the integral (e.g. with the three-point 
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Gaussian quadrature [27j) but results must be equivalent 
after the integration. 

The method used to include the spin-orbit and three- 
nucleon interaction in the propagator, and a detailed de- 
scription of the AFDMC method can be found in Refs. 

MMM- 



C. Importance sampling 

In order to reduce the variance of estimators, the im- 
portance sampling is required. In practice a diffusion 
Monte Carlo calculation is performed using a Green's 
function modified as follows: 



G(R, R', Ar) 



G(R,R',Ar). (33) 



The so called importance function ipi in the above 
equation is often the same as that used for the projection 
of the energy, and is evaluated at the walker configura- 
tion. More precisely we define 



V'/(R) = (^/|R) 



(34) 



In this case the distribution function that is sampled 
in the imaginary time converges to the quantity 



/(R,T^oo) = V/(R)0o(R). 



(35) 



The propagator becomes a shifted Gaussian with a 
modified weight 



1 



'/'j(R) 



-Bi.(R) + Ej.(R') 



Go(R,R',Ar) = 



where D = h? /m is the diffusion constant, and 



Et At 



^l(R) = - 



V^V^7(R) y(R)V^j(R) 
2m V/(R) V^/(R) 



(36) 



(37) 



is the local energy of the system. The additional term in 
the Gaussian is often called drift, so each walker's con- 
figuration is diffused according to 



R' 

where the quantity 



R + DArd 



VV'7(R) 



V'/(R) 



(38) 



(39) 



is the drift term, and 77 is a Gaussian random vector. 

Importance sampling can also be included in the 
Hubbard-Stratonovich transformations that rotate nu- 
cleon spinors. For auxiliary fields importance sampling 
is achieved by "guiding" the rotation given by each 0„ 



operator. More precisely one can consider the following 
identity: 



+ V-A„ATa;„0„ 



+ V-KArXniOn) + v/-A„Atx„ {On ~ (0„)) , 

where the mixed expectation value of the operator (see 
the next subsection for details) is evaluated in the old 
spin configuration: 



{lbl\On\R, S) 



(40) 



This can be implemented by shifting the Gaussian used to 
sample auxiliary fields, and considering the extra terms 
in the weight for branching 



-(a;„-2„)^/2 \/-A„At2;„0„ x„a:„ 



•/2 



where 



X = v^-A„At(0„) . 



(41) 



(42) 



The additional weight term in Eq. |4T]can also be included 
as a local potential, so it becomes 



{HV\R, S) 
{4'i\R,S) 



At 



(43) 



By combining the diffusion, the rotation and all the addi- 
tional factors it is possible to write an explicit propagator 

G(R, R', At) = Go(R, R', At) 



xe 



a2 V^|^j(R,S)| , (^j]V]RS) 
2m \i,j{R,S)\ "I WFWSy 



En ]At 



^i{R,S) \ibiiR',S)\ ' 
where in the choice of the drift term is 
V|^/(i?, S)\ 



\i^i(R,S)\ 



(44) 



(45) 



D. Computation of expectation values 

The projected walker distribution obtained with the 
AFDMC is used to compute expectation values. For a 
generic operator O in the limit t — > 00 the "mixed" ex- 
pectation value is computed as 



)(R)|0|7^t(R)) 
/.o(R)|^t(R)) 



(46) 



We are interested in the expectation value over the 
ground state (/jq. Assuming that i/'t is a good approxima- 
tion of the ground-state a better estimate of the ground- 
state expectation value can be obtained by combining the 
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variational Monte Carlo and the diffusion Monte Carlo 
estimators in this way: 



(0> = 2(0). 



(47) 



where (O)i, is the expectation value computed over the 
variational wave function ipT used as trial wave function. 

The evaluation of the energy of the system is a partic- 
ular case, and can be directly calculated from the pro- 
jected distribution. Since the propagator commutes with 
the Hamiltonian (but this will change in the next section 
when we introduce a constraint), we have 



6o(R)|g|V>T(R)) _ (^T(R)|g|0o(R)) 

(^t(R)|0o(R)) 



,(R)|^t(R)) 



(48) 



The total energy is already the correct value and since it 
does not contain a linear error from the trial function, it 
does not require the extrapolation of Eq. 1471 

The propagator used in our AFDMC calculations is 
written to include only the first eight operators of the 
Argonne interactions. However, in some case, we can 
also evaluate the expectation value of the full Argonne 
Vis Hamiltonian. In light nuclei, the expectation value 
of Argonne Wg is within few percent of Argonne vis ■ It 
is then reasonable to propagate the wave function using 
the Argonne v'g and evaluate the difference between Ar- 
gonne v'g and Vis using the extrapolation of Eg. 1471 This 
procedure was verified in GFMC calculations |13|. and we 
employed this technique in the case of low density where 
Argonne v'g is a very good approximation to Argonne wig. 

More precisely, we evaluate the energy using Argonne 
v'g in the propagator (in addition to the three-nucleon 
interaction), and we add to the total energy the value of 
{vis — v's) evaluated as inU?) We expect this approxima- 
tion to be accurate if this difference is small as in light 
nuclei. 



E. Constrained path and fixed-phase 
approximation 

As described in the above sections, diffusion Monte 
Carlo projects out the ground-state of a given Hamilto- 
nian in terms of the distribution of the walkers. How- 
ever the density of walkers must always be positive 
definite [Slj. For walkers with positive weights, this con- 
dition restricts, in principle, the use of the method to 
that class of problems where the trial wave function is 
always positive or is node-less, such as for a Bose system 
in the ground-state. Algorithms which allow negative 
weights, such as transient estimation [s^ generally have 
exponentially increasing variance. 

One way to deal with Fermionic systems is to set arti- 
ficial boundary conditions between the positive and neg- 
ative regions of the trial wave function. It is possible 
to define a nodal surface where the trial wave function 
is zero and during the diffusion process a walker that 



crosses the nodal surface is dropped; this is the fixed- 
node approximation [s^, [33| . and its application in the 
diffusion Monte Carlo algorithm always gives an upper- 
bound to the true Fermionic ground-state energy. 

In the case of nuclear Hamiltonians or for prob- 
lems where the trial wave function must be complex, 
a constrained-pathjs^, Issl . [SG^ approximation is usually 
applied to avoid the Fermion sign or phase problem. 
The constrained-path method was originally proposed by 
Zhang et al. as a generalization of the fixed-node ap- 
proximation to complex wave functions. In constrained- 
path, walkers are constrained to regions where the real 
part of the overlap with the trial wave function is posi- 
tive. This constrained-path approximation was the orig- 
inal method used to control the phase problem in the 
AFDMC algorithm l^J. More precisely, we have to con- 
sider that even for a complex wave function drift term 
for the coordinates must be real. In the case of the 
constrained-path approximation, a natural choice for the 
drift is 



d 



i?e[^/(R)] 



(49) 



Moreover, in order to eliminate the decay of the signal- 
to-noise ratio it is possible to impose the constrained- 
path approximation, by requiring that the real part of 
the overlap of each walker with the trial wave function 
keeps the same sign. Thus, one can impose 



>0, 



(50) 



where R and R' denote the coordinates of the system 
after and before the diffusion of a time-step. If this con- 
dition is violated, the walker is dropped. This form was 
found to give better results and was employed in the past 
AFDMC calculations [33, [H]. 

An alternative way to control the sign problem is the 
fixed-phase approximation. This method was originally 
proposed by Carlson for nuclear systems [111, and also 
employed for systems whose Hamiltonian contains a mag- 
netic fieldji^. 

We start with the same condition of the reality of the 
drift, and we consider the following expression 



v|V'/(R)l 
l^/(R)l 



(51) 



With this choice the weight for branching becomes 



exp 



v^\^PiCR)\ , yV/(R) 



2m |V'/(R) 



At 



V'/(R) 
|V;7(R)| ^i{R') 

|V'/(R')I V'/(R) 



(52) 



Note that in the above expression there is the usual im- 
portance sampling factor as in Eq. 1331 and an additional 
factor that corrects for the particular choice of the drift. 
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A generic complex wave function can be written as 

^(R) = |V^(R)|e*^(^) , (53) 

where 0(R) is the phase of ip(R.); the factor appearing 
in Eq. [52] can be rewritten as 



IV;, (R) I 

IV/ (R') I V'/(R) 



(R.)] 



(54) 



The fixed-phase approximation constrains the walkers 
to have the same phase as the importance function ipj. 
It can be applied by keeping the real part of the last ex- 
pression. To keep fixed the normalization of the Green's 
function one has an additional factor in the Green's func- 
tion that must be included in the weight: 



(55) 



This can be automatically included by keeping the real 
part of the kinetic energy. In fact: 



Re 



V^IV'/(R)I 
IV'/(R)I 



(56) 



The real part of the kinetic energy includes the additional 
weight term given by the fixed-phase approximation. 

A different derivation to introduce the fixed-phase ap- 
proximation is the following. Let us consider the evolu- 
tion of a complex trial wave function including the im- 
portance sampling: 

■4j}{R)i:{R,T) = J G(R,R',T)V'KR)V'(R',0)dR'. 

(57) 

The quantity ?/'J(R)V'(R, t) is not real and positive defi- 
nite as required, but it is possible to obtain another pos- 
itive density as 

IV'/(R)I ,U(R')-rf,(R)l 



|V'/(R)||V(R,r)| = J G(R,R',r)i^ 



)l 



x|V/(R')IIV'(R',0)MR'. 



(58) 



In this way we impose that the phase of the trial wave 
function is the same of that of ipj. 

Both the constrained-path and the fixed-phase are ap- 
proximations to deal with the Fermion sign problem and 
in principle they should be equivalent if the importance 
function is close to the correct ground-state of the sys- 
tem. 

It is important to note that Carlson et al.i41l| showed 
that within the constrained-path approximation the al- 
gorithm does not necessarily give an upper bound in the 
calculation of energy. This was also observed by Wiringa 
et al. in some nuclear simulations using the GFMC 
technique [12]. It is not guaranteed that our fixed-phase 
calculations give an upperbound. However, in diffusion 
Monte Carlo calculations of the ground-state of quantum 
dots, where both a real or a complex trial wave function 
can be implemented, the fixed-phase approximation gives 
a higher energy that the fixed- node approximation |43j . 



F. Trial wave function 

The trial wave function used as the importance and 
projection function for the AFDMC algorithm has the 
following form: 



V/(R,^) =Fj(R)I?(R,^), 



(59) 



where R = (ri, . . . , tn) represent the spatial and S = 
(si, . . . , sjv) are the spin states of the system. The spin 
assignments Si consist of giving the two-spinor compo- 
nents for each neutron, namely 



a.|T)+&,U), 



(60) 



where ai and bi are complex numbers, and the {It), U)} 
is the neutron-up and neutron-down base. 

The Jastrow correlation function i^j(R) is symmetric 
under the exchange of two particles, and independent of 
spin. Its role is to include the short-range pair correla- 
tions in the trial wave function. The generic form for the 
Jastrow is 



(61) 



where the function /(r) is the solution of a Schrodinger- 
like equation for f(r < d). 



-W'f{r) + v{r)f{r)=\f{r). 
m 



(62) 



where v{r) is the spin- independent part of the nucleon- 
nucleon interaction, the healing distance d < L/2 is a 
variational parameter and L is the size of the box. For 
distances r > d wc impose /(r) — 1. The Jastrow part 
of the trial wave function in the AFDMC case has only 
the role of reducing the overlap of neutrons, therefore 
reducing the energy variance. Since it does not change 
the phase of the wave function, it does not influence the 
computed energy value in projections methods. In all 
the reported results we then fixed d — 2 im or d = L/2 
if L/2 < 2 fm. 

The antisymmetric part of the trial wave function is 
usually given by the ground-state of the non-interacting 
Fermions, which is written as a Slater determinant 



N 



J|0a(r,,s,) 



= Det{(l)a{ri,Si)} , (63) 



where a is the set of quantum numbers of single-particle 
orbitals, and A is the antisymmetrization operator. 

For neutron matter calculations we choose the anti- 
symmetric part as the ground state of the Fermi gas, 
built from a set of plane waves. The infinite uniform sys- 
tem is simulated with N nucleons in a cubic periodic box 
of volume L'^. The momentum vectors in this box are 



2tt 



(64) 
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where a labels the quantum state and rix , Uy and are 
integer numbers describing the state. The single-particle 
orbitals are given by 

(j)airi,Si) = e'''°''''(Xs,m,,a|Si) (65) 



G. Twist-averaged boundary conditions 

Aside from the effect of the phase of the importance 
function employed during the projection in imaginary 
time, the dependence of the energy on the number of 
neutrons is the largest systematic error. Usually one uses 
periodic boundary conditions to reduce finite size effects, 
and simulations are carried out by using a number of neu- 
trons filling closed shells of plane waves. There are still 
sizable errors in the kinetic energy coming from the shell 
structure even at the closed shell filling in momentum 
space (1, 7, 19, 27, 33, 57, ...). In order to estabhsh the 
effect of the finite size of the system due to the kinetic en- 
ergy we imposed twist- averaged boundary conditions [3| 
on the trial wave function. Within periodic boundary 
conditions, the phase, picked up by the wave function as 
a particle makes a circuit across the unit cell, can be cho- 
sen arbitrarily. These more general boundary conditions 
for a wave function are 

i^(ri + Lx, r2, . . .) = e"'-i^(ri, r2, . . .) , (66) 

where L is the side of the simulation cell. The boundary 
condition 9 — gives the usual periodic boundary condi- 
tions, and the more general condition with 9^0, twisted 
boundary conditions. If the twist angle is integrated over, 
the single-particle finite-size effects, arising from shell ef- 
fects in filling the plane wave orbitals, are substantially 
reduced. Integrating over twists averages over the vol- 
ume of k space occupied by the first N Brillouin zones 
of the simulation cell. The occupied region is a convex 
polyhedron that tends to the Fermi surface in the limit 
of infinite system size and has the correct volume at all 
system sizes. The twist averaged kinetic energy must 
approach the exact energy always from above since the 
single-particle kinetic energy is a convex function of k. 

The integration over angles can be achieved in different 
ways, either by modifying the trial wave function during 
the simulation or by performing several simulations using 
different wave functions [441] . In practice, once the density 
of the system is fixed, we consider a grid of different k^- 
vectors 

k„,, = (27rn„ + 6*,)/^ (67) 

within the radius corresponding to the Fermi energy, and 
for each twist angle 9i a simulation is performed. The 
total energy is the given by averaging all the energy ob- 
tained for each wave function. 



H. The Algorithm 

The structure of the AFDMC algorithm consists in the 
following procedures: 

1. Sample the positions and spins, to give \R, S) 
for the initial walkers, from \{'i'j\R, S)\'^ using 
Metropolis Monte Carlo. 

2. Propagate the spatial degrees of freedom as in the 
usual diffusion Monte Carlo with a drifted Gaussian 
for a time step. That is, each walker configuration 
is diffused according to Eq. [38l 

3. For each walker, build and diagonalize the potential 
matrix A'^'^'> . 

4. Loop over the eigenvectors, sampling the cor- 
responding shifted Hubbard-Stratonovich variable 
and update the spinors for a time step. Introduce 
approximate importance sampling of the Hubbard- 
Stratonovich variables, as discussed in the previous 
subsections. 

5. Propagate with the spin-orbit interaction, using 
importance sampling. 

6. Evaluate the real part of the local energy to con- 
strain each walker to have the fixed-phase as de- 
scribed above. This quantity is also stored with 
the corresponding weight to calculate the averaged 
mixed energy. 

7. Iterate from 2 to 6 as long as necessary until con- 
vergence of the energy is reached. 

To evaluate the error bars, block averages are calcu- 
lated and the results combined over different block sizes 
until the blocks become uncorrelated and the error bars 
independent of block size within statistics. 

IV. RESULTS 
A. Test of the fixed-phase approximation 

The AFDMC algorithm combined with the 
constrained-path approximation was previously em- 
ployed by Sarsa et al. to study the neutron matter 
equation of state at zero temperature [1^. In that paper 
the Hamiltonian contained both a realistic Argonne 
v'g two-nucleon and the Urbana IX three-nucleon in- 
teractions; this Hamiltonian is often used to calculate 
properties of both symmetric nuclear matter and pure 
neutron matter. 

The constrained-path AFDMC proved to give very sat- 
isfactory results for neutron matter calculations with a 
two- and three-nucleon interactions, but some problems 
were encountered in the evaluation of the spin-orbit con- 
tribution. The inclusion of spin-backflow correlations 
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reduced the discrepancies. A detailed study consider- 
ing a pure nucleon-nucleon interaction emphasized the 
problem of constrained-path AFDMC in dealing with 
the spin-orbit interaction ^455. A similar behavior was 
found by comparing the constrained-path AFDMC with 
the GFMC evaluation for the energy of 14 neutrons in a 
periodic boxQ. When using the same Hamiltonian with 
the same box truncation used in GFMC calculations of 
Ref. , the constrained-path AFDMC overestimated the 
energy of 14 neutrons with an Argonne v'g interaction. 

The AFDMC with the fixed-phase approximation over- 
comes the discrepancies previously observed in the es- 
timates of the spin-orbit contribution to the total en- 
ergy, as shown in tabic |T1 Without tail corrections 
the constrained-path AFDMC energy of 14 neutrons at 
p =0.16 fm-3 is 20.32(6) MeV compared to 17.00(27) 
MeV given by unconstrained GFMC[^, while the fixed- 
phase AFDMC energy is 17.67(5) MeV, within 3%, and 
in much better agreement with unconstrained GFMC. 



p [fm-3] 


FP-AFDMC CP-AFDMC CP-GFMC UC-GFMC 


0.04 
0.08 
0.16 
0.24 


6.75(7) 6.43(01) 6.32(03) 
10.29(1) 10.02(02) 9.591(06) 
17.67(5) 20.32(6) 18.54(04) 17.00(27) 
27.7(5) 30.04(04) 28.35(50) 



TABLE I: Fixed phase (FP-AFDMC) energies per particle 
of 14 neutrons interacting with the Argonne v'g interaction 
in a periodic box without the inclusion of finite size efi'ects 
at various densities. The constrained-path (CP-AFDMC) of 
Ref. [1^, the constrained-path (CP-GFMC) and the uncon- 
strained (UC-GFMC) GFMC of Ref. 9] are also reported for 
a comparison. All the energies are expressed in MeV. 

For higher densities reported in table |T] it should be 
noted that the constrained-path GFMC significantly dif- 
fers from unconstrained GFMC because the Fermion sign 
problem becomes more severe and the unconstrained en- 
ergy estimation has larger fluctuations. The convergence 
can be hard to reach because the imaginary time evolu- 
tion of the energy can be carried out only for very small 
steps. These reasons could introduce some spurious ef- 
fects limiting the accuracy of GFMC for the neutron mat- 
ter calculation to densities below 0.08 fm~^ 0. 

Preliminary results for the ground-state calculation 
of neutron drops by means of the fixed-phase AFDMC 
show that the spin-orbit contribution is now in agree- 
ment with the GFMC results [S^l- Using the same Hamil- 
tonian, previous constrained-path AFDMC calculations 
predicted a spin-orbit splitting (SOS) in "^n neutron-drop 
about a half the GFMC result [11]. Instead the fixed- 
phase AFDMC estimate is in excellent agreement with 
the GFMC one, 30 1 , also for neutron drops containing up 
to 13 neutrons |46|. 

The improvement yielded by using the fixed-phase 
approximation rather than the constrained-path is also 
evident in the comparison of the fixed-phase AFDMC 
with the available constrained-path AFDMC using spin- 



backflow correlations. In table [IT] we report all the avail- 
able calculations computed within the constrained-path 
approximation compared to the fixed-phase one. The 
corrections included are only due to the truncation of the 
nucleon-nucleon interaction as in the old calculations [28|. 



p [im ^] 


FP-AFDMC(14) CP-AFDMC(14) JSB-AFDMC(14) 


0.12 


14.52(5) 


14.80(9) 


0.16 


19.03(7) 


19.76(6) 


0.20 


24.49(5) 


25.23(8) 


0.32 


46.60(8) 


48.4(1) 46.8(1) 


p [fm-»] 


FP-AFDMC(66) CP-AFDMC(66) JSB-AFDMC(66) 


0.12 


15.04(8) 


15.26(5) 


0.16 


20.14(5) 


20.23(9) 


0.20 


26.21(5) 


27.1(1) 


0.32 


52.47(4) 


54.4(6) 52.9(2) 



TABLE 11: Fixed phase (FP-AFDMC) energies per particle of 
14 and 66 neutrons interacting with the Argonne tig -|- Urbana- 
IX interaction in a periodic box at various densities com- 
pared with the available constrained-path (CP-AFDMC) ones 
of Ref. 111]. The constrained-path AFDMC results us- 
ing a Jastrow-Slater-backflow (JSB) wavefunction 47] are also 
shown. In order to make the comparison possible, the finite 
size effect due to the truncation of nucleon-nucleon interac- 
tion was included, while that of Urbana-IX was omitted. All 
the energies are expressed in MeV. 



B. Equation of state of neutron matter 

We employed the fixed-phase AFDMC method to 
study neutron matter by simulating different numbers 
of neutrons interacting with the Argonne v'g potential, 
including finite-size corrections as described in Ref. [281] . 
All of the fixed-phase AFDMC results are reported in ta- 
ble mil which shows the energy per neutron of neutron 
matter for different densities by varying the number of 
neutrons. 

Some finite-size effects are present, as can be deduced 
by observing the energies for different numbers of neu- 
trons. The same behavior is followed at each density, 
and £;(38) < E{U) < E{66). This trend directly follows 
the kinetic energy oscillations of N free Fermions which 
for iV = 38 is lower than either = 14 or iV = 66. 

For neutron matter the Urbana-IX three-nucleon force 
reduces to a pairwise spin interaction modulated by the 
spectator neutron as explained in Refs. f2^,'3^ which can 
be easily included into the propagator. Finite-size correc- 
tions due to the Urbana-IX can be included in the same 
way as for the nucleon-nucleon interaction, although their 
contribution is very small compared to the potential en- 
ergy. Their effect is appreciable only for a small number 
of particles and at large density, i.e. if the size of the 
simulation box is small. 

All the fixed-phase AFDMC results of 14, 38, 66 and 
114 neutrons interacting with the Argonne v'g Urbana-IX 
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p [fm-«] 


E/N (14) E/N (38) E/N (66) 


0.12 


12.08(5) 


11.18(4) 


12.65(4) 


0.16 


14.87(9) 


13.50(5) 


15.43(3) 


0.20 


17.6(1) 


16.10(4) 


18.27(5) 


0.24 






21.56(5) 


0.28 






25.05(6) 


0.32 


27.2(1) 


25.2(1) 


28.93(7) 


0.36 






33.05(6) 


0.40 






37.15(8) 


0.48 






46.7(1) 


0.56 






57.64(9) 


0.64 






69.90(8) 


0.80 


91.5(2) 


89.2(2) 


97.4(1) 



TABLE III: Fixed phase AFDMC energies per particle of 14, 
38 and 66 neutrons interacting by the Argonne potential 
in a periodic box at various densities. The finite-size effects 
due to the nucleon-nucleon truncation are included. All the 
energies are expressed in MeV. 



p [fm-^] 


E/N(14) 


E/N(38) E/N (66) E/N (114) 


0.12 


14.77(7) 


13.68(3) 


15.18(2) 


16.05(4) 


0.16 


19.41(7) 


18.32(4) 


20.04(2) 


21.31(4) 


0.20 


25.05(7) 


24.06(4) 


26.13(4) 


27.82(5) 


0.24 


31.74(6) 




33.64(4) 




0.28 


39.79(3) 




42.51(3) 




0.32 


48.61(5) 


48.76(6) 


51.84(2) 


55.13(6) 


0.36 


60.03(5) 




64.89(5) 




0.40 


72.38(5) 




78.59(6) 




0.48 


102.74(5) 




111.69(9) 




0.56 


139.8(1) 




152.81(2) 




0.64 






202.19(9) 




0.80 




320.3(1) 


328.19(6) 





TABLE IV: Fixed phase AFDMC energies per particle of 
14, 38, 66 and 114 neutrons interacting by the Argonne 
iig+Urbana-IX potential in a periodic box at various densi- 
ties. Note the difference with the values of table |ll] due to 
the different treatment of finite size effects, that, in this case, 
include two- and three-nucleon interaction contributions. All 
the energies are expressed in MeV. 

Hamiltonian, and including all the finite size effects due 
to the truncation of two- and three-nucleon interactions 
are summarized in table Hvl Important finite size effects 
are still present. The value closest to the thermodynamic 
limit is that for 66 neutrons, because the free Fermi gas 
energy of this particular system is very similar to that of 
the infinite one. However the difference of the energy of 
66 and 114 neutrons is always within 6-7%. This behav- 
ior was also observed in a study of finite-size effects using 
the periodic box FHNC technique [4^], and the analysis 
of Sarsa et ahfl^ suggests that the energy of the system 
in the infinite limit is somewhere between the energies of 



FIG. 1: (color online) Convergence of the computed energy 
at p=0.32 fm"'' as a function of neutrons in a box within 
the grid twist- averaging method (TABC) described in the 
text with ten twists: the Argonne iig-|-Urbana-IX Hamilto- 
nian were considered. The equation of state is compared with 
the fixed-phase AFDMC calculations with periodic boundary 
conditions (PBC) shown by solid lines. 



66 and 114 neutrons. 

In order to better understand the finite size effects due 
to the kinetic energy, we repeated several simulations by 
imposing twist-averaged boundary conditions in the trial 
wave function. The results are displayed in Fig. [1] where 
we reported the energy obtained by averaging all the re- 
sults using sets of ten twist angles in each dimension. The 
different behavior of the energy as a function of the num- 
ber of neutrons using periodic or twist-averaged bound- 
ary conditions is well evident. As expected the effect 
of twist averaging is to reduce the jumps of the energy 
as a function of TV given by periodic boundary condi- 
tions. Then the extrapolation to the infinite limit of N 
is better evident using twist averaging. However it is re- 
markable that the energy of 66 neutrons computed using 
either twist averaging or periodic boundary conditions is 
almost the same. This essentially follows the fact that 
the kinetic energy of 66 fermions approaches the infinite 
limit very well. In addition, the twist averaging could 
be very useful to simulate systems for which an arbitrary 
number of Fermions is needed 49]. 

In Fig. [2] we plot the fixed-phase AFDMC equation 
of state, obtained with the energy of 66 neutrons, and 
the calculation of Akmal et al. of Ref. 50], where the 
Argonne fis interaction combined with the Urbana-IX 
three-nucleon interaction was considered. As it can be 
seen both the Argonne and wig give an equation of 
state showing essentially the same behavior, with a dif- 
ference in the energy that is similar throughout the con- 
sidered range of densities. The addition of the three- 
nucleon interaction increases the differences between the 
AFDMC and the Akmal et al., in particular at higher 
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P [fm '] 



FIG. 2: (color online) The fixed-phase AFDMC equation 
of state evaluated by simulating 66 neutrons in a peri- 
odic box; the Argonne Wg (AV8')and Argonne iig-l-Urbana- 
IX (AV8'-|-UIX) Hamiltonians were considered. The equa- 
tions of state are compared with the variational calculations 
of Ref. [5^ using the Argonne vis (AV18) and the Argonne 
^18 -l-Urbana-IX (AV18-UIX) Hamiltonians. See the legend 
for details. 



The AFDMC results have been fitted with the follow- 
ing functional form: 

^{p) = apP + cp\ (68) 

where E/N is the energy per neutron in MeV as a func- 
tion of the density in fm~^. The parameters of the fit 
for both Argonne and the full Argonne Wg-|-Urbana-IX 
Hamiltonian are reported in table |Vl We also tried to 
use the functional form of Ref [5lj where /3 = 1. We had 
a worse but the equation of state and the pressure as 
a function of the density does not change in a significant 
way. 



Hamiltonian 


a c /3 7 


AV8' 
AV8'-fUIX 


23.0 115.7 0.37 1.87 
32.6 507.8 0.48 2.375 



TABLE V: The parameters of Eq. |68l fitting the equa- 
tion of state computed with the full Argonne «g-|-Urbana-IX 
(AV8'+UIX) Hamiltonian and with the nucleon-nucleon in- 
teraction only (AV8'). The parameters a and c are expressed 
in [MeV/fm-^]. 



densities, implying a strong difference in pressure and 
compressibility. 

The Argonne Wg interaction should be more attractive 
than Argonne wig as shown in light nuclei and in neu- 
tron drop calculation flS*!. The result shown in Fig. ^ 
where we compare Argonne v'g result with Akmal et al.'s 
Vis values, do not show this. We believe this is indica- 
tive of systematic errors in the FHNC/SOC calculations. 
The fixed-phase AFDMC has proved to be in very good 
agreement with the CFMC results for light nuclei [ll|, 
and also with the CFMC results for 14 neutrons. On the 
other hand the fixed-phase AFDMC calculation of the 
nuclear matter suggested that the FHNC/SOC approxi- 
mation could miss important contributions, in particular 
those coming from the neglected elementary diagrams in 
the FHNC summation il2|. In the case of Akmal et al. 
calculations the energy is computed by means of a clus- 
ter expansion for which no evidence of convergence can 
be provided. The addition of the Urbana-IX three-body 
interaction to the Hamiltonian increases the differences 
between the AFDMC results and of Akmal et al. ones, 
and, again, this confirms that the variational technique 
based on the cluster expansion gives a lower energy be- 
cause it neglects important contributions. However, we 
stress the fact that in the case of neutron matter the con- 
tribution of the tensor-r force is small compared to the 
other channels of the interaction. For such reason the 
calculation of the energy within traditional variational 
techniques based on FHNC/SOC or cluster expansion 
could be more accurate for pure neutron matter with- 
out protons. This is not true in dealing with nuclear 
matter where the effect of tensor-r is most important, as 
confirmed in ref. [l2|. 



C. Argonne ng and vis interactions 

As described in the above sections, in most cases the 
Argonne vis result is evaluated as a perturbation of the 
Argonne VgilS]. The assumption is reasonable since the 
Argonne v'g potential contains most of the contributions 
of vig, potential and was obtained with a reprojection by 
keeping only the most important terms. However the 
operators appearing in Argonne vis and not in Argonne 
Dg are not exactly included in the CFMC calculations. 
The imaginary-time CFMC evolution is performed using 
Argonne v'g and the energy is calculated perturbatively 
in the difference between Wg and wig, which for nuclei is 
a fraction of an MeV. 

This method is also used in the FHNC/SOC calcu- 
lation where only the lowest order two-body nucleon- 
nucleon correlations are included in the variational wave 
function [50|. However there is no reason to believe that 
such a calculation gives an upper bound to the true en- 
ergy, and this approximation may not be good, particu- 
larly for higher densities. 

When using a propagator including the Argonne Wg 
potential the difference between the energies computed 
using Argonne Ug and uig is actually very small. For in- 
stance, for 14 neutrons at p <0.12 fm""^ the difference be- 
tween Argonne Wg and wig is less than 2 MeV per neutron 
and is 2.7 MeV and 5.1 MeV for densities of 0.16 fm"^ 
and 0.20 fm~^ respectively. On the other hand, a plain 
truncation of Argonne vis in the propagator leads to huge 
energy differences in the two estimates. At p=0.12 fm~'^ 
the energy of 14 neutrons with the Argonne Ug -I- Urbana- 
IX Hamiltonian is 14.12 MeV, while it is 3.60 MeV for 
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Argonne wig+Urbana-IX. This means that the extra vis 
terms cannot be thought of as a small correction to Ar- 
gonne Ug, at least in this range of densities. However, 
at p <0.04 fm~^ the difference between Argonne Wg and 
Vis is a few percent of the total energy, so we can safely 
evaluate this difference as a perturbation using the Wg 
propagator. 

In the very low-density regime the neutron matter is a 
superfluid gas, and a trial wave function written in a BCS 
form including explicitly the pairing between neutrons is 
needed[52, 53]. However, we expect that the expectation 
value of the Argonne Wg interaction to be of order of that 
of Argonne uig both in the superfluid and in the normal 
phase. Here we are only interested to a qualitative study 
of the difference between Argonne Ug and wig, thus a wave 
function as presented in section IIII Fl was used, rather 
then that of Ref. 

It is interesting to focus on the equation of state of 
neutron matter in the low-density regime and in the nor- 
mal phase. The Argonne wig-|-Urbana-IX Hamiltonian 
as described was used. The range of p <0.04 fm~^ is 
particularly relevant in the study of properties of the in- 
ner crust of neutron stars. The very low density neu- 
tron matter approaches a regime which is almost univer- 
sal, and is analogous to, for instance, cold atoms0|. In 
this regime many-body techniques can be compared and 
calibrated [H, [Hgl . 





E/N 


(AV8'-AV18) 


3.377 X 10"^ 


0.089(1) 


-0.00197574 


2.702 X 10"'' 


0.367(2) 


0.0002776 


0.002162 


1.289(2) 


0.002525 


0.007295 


2.606(4) 


0.021712 


0.01729 


4.277(7) 


0.082534 


0.03377 


6.197(2) 


0.30802 



TABLE VI: Fixed phase AFDMC energies per particle 66 
neutrons interacting with the Argonne Ug -fUrbana-IX inter- 
action in a periodic box at various densities. The difference 
between the t;g and the vis interactions (AV8'-AV18), evalu- 
ated with an extrapolation as described in Sec. IIII Dl is also 
reported. All the energies are expressed in MeV. 

We report the energy of 66 neutrons in a periodic 
box in table IVII The Hamiltonian uses the Argonne 
z)g-|-Urbana-IX potential, and the potential is corrected 
for finite size effects. The difference between the Argonne 
v's and the fig interactions was extrapolated as described 
in Sec. IIII Dl As it can be seen, at such densities the Ar- 
gonne Wig contribution is similar to Ug , of the order of few 
percent with respect to the total energy. We believe that 
such a difference would be even smaller if the full Argonne 
Wig was implemented in the propagator, and then eval- 
uated without the extrapolation described in the above 
sections. 

The equation of state of neutron matter in the low- 
density regime is reported in Fig. [31 where the energy 
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FIG. 3: (color online) The equation of state of neutron mat- 
ter in the low-density regime. The Argonne Ug (AV8') and 
v\s (AV18) interactions combined with the Urbana-IX (UIX) 
three-nucleon interaction were considered as indicated in the 
legend. 



per neutron as a function of the density is calculated 
both with the AV8' and with the AV18 nucleon-nucleon 
interaction combined with the Urbana-IX three-nucleon 
interaction. As it can be seen, the difference between the 
two Hamiltonians considered is very small in this regime. 
The Argonne Wg and wig combined with the Urbana-IX 
essentially give the same energy, and only small devia- 
tions are present when the density increase above wO.OlS 
fm~'^. This result is confirmed by the fact that in such 
a regime the neutron-neutron interaction is dominated 
by the S channel that in the Argonne Wg is the same 
as in uigfi^]- There is a small trend that the energies 
are sensibly higher at p <0.015 fm~'^. Other many-body 
calculations are in general not in agreement and present 
very different behaviors in this regime [58i|. 



V. CONCLUSIONS 

We accurately calculated the equation of state of neu- 
tron matter using the auxiliary field diffusion Monte 
Carlo method. We started from a non-relativistic nu- 
clear Hamiltonian containing two- and three-nucleon po- 
tentials. The AFDMC algorithm suffers from the usual 
fermion sign problem present in all fermionic Monte 
Carlo calculations, and we find that the fixed-phase 
approximation used to control it to be more accurate 
than the previously used constrained-path approxima- 
tion. In particular, in this work we demonstrated that 
the fixed-phase AFDMC overcomes the problems encoun- 
tered when dealing with the spin-orbit interaction. The 
Urbana-IX three-body force is included in the fixed-phase 
AFDMC calculation without any perturbative evalua- 
tion because it is naturally included in the Green's func- 
tion used for the propagation. The fixed-phase AFDMC 
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reveals some problems of the variational cluster sum- 
mation (or FHNC/SOC) technique just highlighted in 
the nuclear matter calculation with a simple ug-like 
interactionfl^. 

We computed the equation of state of neutron mat- 
ter using a modern, but still simplified, nucleon-nuclcon 
interaction combined with a realistic three-nucleon inter- 
action in the regime of interest for predicting the proper- 
ties of neutron stars, and we found some deviations with 
respect to past variational calculations based on cluster 
expansion, in particular at high densities. The differ- 
ence between the Argonne Wg-|-Urbana-IX Hamiltonian 
and that containing the more sophisticated Argonne vis 
interaction was perturbatively evaluated in the low den- 
sity regime, where the equation of state is useful to con- 
strain properties of the inner crust of neutron stars. Our 
equation of state can also be useful to compare the wide 
range of Skyrme forces used to study the neutron matter. 

We are working to include the full Argonne vis inter- 



action in the two-body part of the Hamiltonian, and we 
are investigating the effect of the more complex Illinois 
three-nucleon forces. The effect of those forces in neutron 
drops and in neutron matter will be a subject of future 
work. 
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